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One of the most striking features of novel 2D semiconductors (e.g., transition metal dichalcogenide 
monolayers or phosphorene) is a strong Coulomb interaction between charge carriers resulting in 
large excitonic effects. In particular, this leads to the formation of multi-carrier bound states upon 
photoexcitation (e.g., excitons, trions and biexcitons), which could remain stable at near-room 
temperatures and contribute significantly to optical properties of such materials. In the present 
work we have used the Path Integral Monte Carlo methodology to numerically study properties of 
multi-carrier bound states in 2D semiconductors. Specifically, we have accurately investigated and 
tabulated the dependence of single exciton, trion and biexciton binding energies on the strength 
of dielectric screening, including the limiting cases of very strong and very weak screening. The 
results of this work are potentially useful in the analysis of experimental data and benchmarking of 
theoretical and computational models. 


I. INTRODUCTION 

Last decade has witnessed an ever-growing interest of 
the scientific and engineering communities in atomically 
thin two-dimensional (2D) materials [IMj. This is due 
to both the unique properties of such materials from the 
basic science perspective 01, as well as their tremen¬ 
dous technological potential HIM]- While graphene - 
a flagship of this class of materials - has attracted most 
of the attention to date, it lacks an electronic bandgap 
and is thus unsuitable “as is” for multiple applications 
including photovoltaics and electronics 1 [TO]. The re¬ 
alization of this fact has led to significant efforts spent 
during the last few years on searching for “graphene with 
a bandgap”, i.e., atomically thin semiconductor materi¬ 
als. Of many materials identified as potential semicon¬ 
ductor alternatives to graphene, one of the most promis¬ 
ing classes of materials is monolayer transition metal 
dichalcogenides (ML-TMDC) [TTJ [T2] . Other possible 
candidates include phosphorene H3HH , stanene mum 
and ultrathin organic-inorganic perovskite crystals nu. 

A prominent feature of 2D semiconductors (2DS) is 
the strong Coulomb interaction between charge carri¬ 
ers, as compared to their parent 3D materials (e.g., bulk 
TMDC). This effect arises from the so called dielectric 
confinement - the dielectric screening of the Coulomb 
interaction between two point charges within a 2D ma¬ 
terial becomes weak at distances larger than some effec¬ 
tive thickness of the material | [20l423] . In other words, 
the strength of dielectric screening becomes distance- 
dependent, i.e., spatially non-local. For example, this 
effect is pronounced in ML-TMDC, leading to the forma¬ 
tion of tightly-bound excitons with binding energies up to 
~ 0.4 —0.7eV i2T| - 126| . This is very different from conven¬ 
tional 3D semiconductors, where the optical response at 
room temperature is dominated by independent charge 
carriers m- Such large exciton binding energies in 2DS 
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suggest the possibility of existence of other stable bound 
states consisting of a larger number of charge carriers 
(e.g., trions, biexcitons) at near-room temperatures [281 . 
Indeed, trions - bound states of a single exciton and an 
extra charge carrier - with binding energies on the order 
of ~ 20 — 45 meV in ML-TMDC and up to ~ 190 meV in 
few-layer phosphorene have been recently observed PSH- 
[33] • Very recent evidence of stable bound states of two 
electrons and two holes - biexcitons - with binding ener¬ 
gies of ~ 45 — 70 meV in ML-TMDC has been reported 
in Refs. [31M33| . 

The term “excitonic effects” has been traditionally ap¬ 
plied to features of the optical properties of semicon¬ 
ductors, which go beyond the simple picture of non¬ 
interacting charge carriers. In particular, the Coulomb 
interaction between carriers leads to the formation of 
multi-carrier bound states out of free charge carriers. 
Theoretically, excitonic effects have mostly been studied 
at the level of a single exciton using electronic structure 
theory tools BUB1[2S1|37H1I]. Experimentally observed 
large exciton binding energies were largely reproduced 
and the detailed structure of excitonic states was studied. 
In particular, it has been demonstrated that even without 
any external dielectric screening (e.g., due to a solvent or 
substrate) the size of an exciton in ML-TMDC is much 
larger than the size of a single unit cell donum]- Under 
these conditions, quasimomenta of interacting electrons 
and holes are distributed very close to the positions of 
the respective band edge extrema in the Brillouin zone, 
see for example the inset of Fig. 3(b) in Ref. [38] • This 
observation justifies the use of the effective mass approx¬ 
imation in the analysis of excitonic effects in ML-TMDC 

B2HM1I3S]- 

Despite enormous interest, however, theoretical and 
computational studies of multi-carrier bound states be¬ 
yond a single exciton (e.g., trions, biexcitons) in 2DS are 
rather sparse to date. To the best of our knowledge, 
the first theoretical analysis of single excitons and trions 
in novel 2DS explicitly taking into account the spatially 
non-local character of dielectric screening was performed 
in Ref. [23], followed by Refs. [22] (42j (43|. In partic¬ 
ular, trion binding energies were found by means of a 
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simple few-parameter variational wave function [23J I43j . 
Trion excitations in these systems have also been studied 
by means of the fractional dimensional space approach 
[44| . and the time-dependent density-matrix functional 
theory |45| . but no non-local screening effects were ex¬ 
plicitly taken into account in these works. Finally, an 
accurate analysis of single exciton and trion binding en¬ 
ergies in the limit of a very strong non-local 2D screening 
has been reported very recently |46j . 

The goal of the present work is to accurately study the 
excitonic effects in 2DS at the arbitrary strength of the 
non-local 2D screening. Within the effective mass ap¬ 
proximation, we demonstrate how a general problem of 
finding the energy of a bound state of a specific number 
of charge carriers can be reduced to finding an unknown 
function of two parameters: the electron-to-hole mass ra¬ 
tio and the screening length. We then evaluate this func¬ 
tion numerically using the Path Integral Monte Carlo 
(PIMC) methodology )47j . This specific flavor of the 
broad class of quantum Monte Carlo methods is among 
the most general methods to study interacting quantum 
many-body systems. It does not require the assump¬ 
tion of positivity of a wave function (diffusion quantum 
Monte Carlo) @5], or choosing a wave function ansatz 
(variational quantum Monte Carlo) {49], and has been 
successfully applied to analyze the electronic structure 
of low-dimensional semiconductors |50H53) . In particu¬ 
lar, binding energies of excitons, trions and biexcitons in 
GaAs/AlGaAs quantum wells have been accurately eval¬ 
uated Ell- 

Using PIMC, we evaluate the exciton, trion and biex- 
citon binding energies for a range of values of the screen¬ 
ing length assuming identical isotropic electron and hole 
effective masses. This constitutes the main result of the 
paper, since a straightforward interpolation and rescaling 
of the obtained numerical results (compiled in Table [I]) 
can be used to evaluate binding energies of single exci¬ 
tons, trions and biexcitons in arbitrary 2DS with not too 
different isotropic effective masses (e.g., ML-TMDC or 
stanene). To the best of our knowledge, this work is the 
first to numerically exactly treat trions and biexcitons 
in 2DS, albeit within the effective mass approximation. 
Numerical results, obtained here, can be directly used to 
analyze experimental data and benchmark other theoret¬ 
ical and computational approaches to excitonic effects in 
2DS. 

The paper is organized as follows. The general Hamil¬ 
tonian of the system of interacting charge carriers is in¬ 
troduced in Sec. |TT] The dependence of the energy of a 
multi-carrier bound state on the electron to hole effective 
mass ratio and the strength of 2D screening is discussed 
in Sec. m Brief description of the PIMC technique and 
its application to finding energies of multi-carrier bound 
states is given in Sec. m Section |V| contains the discus¬ 
sion of the numerical results. An example of how these 
numerical results could be applied in a realistic situation 
is given in Section m Section m contains the conclu¬ 
sion. 


II. GENERAL THEORY 


Within the effective mass approximation, the Hamilto¬ 
nian for an arbitrary number of interacting charge carri¬ 
ers (electrons and holes) inside a 2D semiconductor (2DS) 
material can be written as [201 1231 [54] 



where pf 



dvl 


and TOj are the isotropic ef¬ 


fective masses. We assume only two types of charge car¬ 
riers with effective masses m e (electrons) and rrih (holes) 
and respective charges q e = —e and qh = e, where e = |e| 
is the magnitude of the electron charge. Number of carri¬ 
ers of each type is encoded by S. In this work, we restrict 
ourselves to 6 = X, T+, T_ or XX, which corresponds to 
a single exciton (e + h), positive trion (e + 2 h), negative 
trion (2e + h), or biexciton (2e + 2h), respectively. Ac¬ 


cordingly, encodes 


a summation only over carriers 


specified by S. 

In this work, the energies of single carriers are defined 
relative to corresponding band edges, so that the 2DS 
bandgap does not explicitly enter Eq. ([I]). Consequently, 
the obtained energies of multi-carrier bound states are 
negative, and, for example, the energy of a single exciton 
obtained from Eq. ([!]) is not the energy of an exciton peak 
directly observed in a photoluminescence experiment. In¬ 
stead, such a negative energy corresponds to an energy 
release upon the formation of a multi-carrier bound state 
out of free charge carriers. 

The non-locally screened Coulomb potential in Eq. 0 
is given by 


V(x) = ^[H 0 (x)-Y 0 {x)], ( 2 ) 


where H 0 (x) and Yq(x) are the Struve function and the 
Bessel function of the second kind, respectively [55]. The 
effective dielectric constant of the environment is defined 
as k = (ki + K 2 )/ 2 , where and «2 are the dielectric 
constants of two materials that the 2DS is surrounded 
by. Screening length is denoted by r^. In the case of a 
thin semiconductor layer of a finite thickness d and an 
isotropic dielectric constant kq, the screening length is 
evaluated as r 0 = dn 0 /2n [2D] ■ In the case of a truly two- 
dimensional material (e.g., atomically thin), the screen¬ 
ing length is given by [54] 


A) = 2ttx2 d/k, (3) 

where Xi d is the 2D polarizability of the material. 
Asymptotic expansions of V(x) at large and small x are 

V(x) » 1/x (4) 

and 

V(x) ss log(2/x) - 7 , (5) 
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respectively, where 7 = 0.57721... is Euler’s constant 
m l54| . The former limiting case corresponds to the 
Coulomb interaction unaffected by the dielectric polar¬ 
ization of 2DS as most of the electric field lines between 
two distant charges go outside of the 2DS. In the opposite 
limit, Xij <C ro, most of the electric field lines are con¬ 
fined within the 2DS. In this case, the Coulomb potential 
becomes logarithmic. 

The dimensional analysis of Hamiltonian ([lj with only 
a single type of carrier for each charge sign results in the 
following universal expression for the ground state energy 

Eg = -2E 0 £s ( a , f 0 ), (6) 

where £g is some yet unknown dimensionless function of 
two continuous real arguments, that does also depend on 
the number of electrons and holes in a system via the 
parameter 5 = X, T+, T_ and XX. The arguments are 
a = m e /rrih - ratio of electron and hole effective masses, 
and ro = ro/a - screening length normalized by the Bohr 

t2 

radius. The Bohr radius is introduced as a = 7, where 

/i = (m ” 1 + m )) 1 ) -1 is the reduced mass. The Bohr 
energy is E 0 = 

III. GENERAL PROPERTIES OF £ s (cr, f 0 ) 

General properties and the asymptotic behavior of 
£g(a,r 0 ) as a function of mass ratio a and the dimen¬ 
sionless 2D screening length f 0 can be understood with¬ 
out performing numerical calculations. 


which can be directly obtained by differentiation of both 
sides of Eq. ([7|. For trions, this equation means that if, 
for example, the energy of the positive trion decreases 
with er at fixed /r and a ~ 1 , then the energy of the nega¬ 
tive trion necessarily increases with the same magnitude 
of the derivative. For charge-neutral multi-carrier states, 
i.e., excitons and biexcitons in this work, we have S = S, 
and, therefore 


d_ 

da 


£s{cr,r 0 ) 


= 0 . 

CT— 1 


(9) 


Since the energy of a single exciton does not depend on 
the electron or hole effective mass independently, but 
only through the reduced mass /x, £x(a,fo) does not 
actually depend on a, so that Eq. § is redundant for 
S = X. However, £g(a, r 0 ) does in general depend on the 
mass ratio for bound states of more than two carriers. 
Consequently, Eq. (|9| suggests a somewhat non-trivial 
result for biexcitons - weak dependence of £xx(c, fo) on 
a when electron and hole effective masses are similar. 

In the opposite limit, i.e., when the mass ratio a be¬ 
comes very small, one has m e ~ n and mh « \i/a m e . 
Under these conditions, the Born-Oppenheimer approxi¬ 
mation can be applied by treating holes as classical par¬ 
ticles. Then, the first quantum-mechanical correction to 
the relative motion of holes would be to treat that mo¬ 
tion as vibrational normal modes of an “electron-hole 
molecule” with frequencies ui ~ m h . The zero-point 
energy of such vibrational normal modes leads to a non- 
analytic asymptotic behavior of £g(a,ro) at a —> 0 


A. Dependence on a 


£ 5 ( 0 , f 0 ) - £g(a,r 0 ) oc a 1/2 . (10) 


First of all, it is important to note that since the Hamil¬ 
tonian 0 does possess the charge conjugation symmetry 
(C-symmetry), it is sufficient to specify £s{a,fo) only in 
the range of a € (0,1). Indeed, due to the C-symmetry 
one always has 


£s(cr,r 0 ) = £g(l/a,r 0 ), (7) 

where 6 specifies a multi-carrier state obtained from state 
S by charge conjugation. For example, if <5 = T+, then 
5 = X_. Therefore, knowing the energy of T + at a < 1 
immediately gives the energy of T_ at cr > 1 and vice 
versa. Substitution cr —► 1/a naturally arises from the 
charge conjugation of cr = m e /nih. Two important con¬ 
sequences of Eq. |7| are as follows. The first rather triv¬ 
ial one is that energies of positive and negative trions 
are the same at a = 1. Henceforth, we will use desig¬ 
nation S = T for trions when electron and hole masses 
are identical and therefore the overall charge of a trion is 
irrelevant for energetics. 

The second consequence is that 


d_ 

for 


Zsiv, fo) 


d_ 

for 


£a(c,f 0) 


( 8 ) 


These considerations, however, require the presence of 
more than a single hole. Therefore, we expect Eq. (10) 
to be accurate for T+ and XX in this work. For XL, we 
expect a more regular behavior at cr —> 0 . 


B. Dependence on r0 


In this subsection, we consider the behavior of £g(a, ro) 
in the limit of very weak (fo —> 0 ) and very strong 
(f 0 —¥ 00 ) 2D screening. In the former limit, the screen¬ 
ing becomes strictly local as it originates only from the 
3D environment (via k) and not from 2DS. In this limit, 
the exact analytical solution of the single exciton prob¬ 
lem exists yielding £x(a,r o)| f _^ 0 = 1 El|. Hence the 
choice of the prefactor of 2 in Eq. (| 6 |. In the case of 
a trion or biexciton, the exact analytical expression for 
£g(a,0) is not available, but multiple approximate ana¬ 
lytical, as well as approximate numerical and numerically 
exact results have been obtained to date [HI 311 f5M55] . 

In the opposite limit of a very strong 2D screening, 
fo —> 00 , the Coulomb potential becomes logarithmic 


a—1 


a—1 
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[see Eq. ([5])], so that Hamiltonian (|T]) can be rewritten as 


H s = — 

tia 


1 v— > 7 ' 1 * 

— A 2 . 

9 * ^ m • 


m. C log 2 A _ 7 

r o V Xij 

i,j<i 


(ii) 

where all the variables inside the square brackets are 
made dimensionless: rhi = rrii/n , Xi = Xi/a, qi = 
qi/e = ±1. We further rescale coordinates by introducing 
= oti/r g/ 2 , which yields 



H s = — 

KT 0 

hg + Eg , 

(12) 

where 

1 1 * 

hs = ~2^ rh^ 

- 'Mjiog 

(13) 


i 

i,j<i 


and 




E! = 

e 2 N e + N h - {N e 

-Nh) 2 (1 , 

- ( -log4r 0 - 

7 )’ 

kxq 2 


(14) 

where N e and Nh are the numbers of electrons and holes, 
respectively, within a multi-carrier bound state. Since 
the effective masses in atomic units can be expressed as 
m e = 1 + a and fhh = 1 + ct - 1 , Hamiltonian in Eq. (13) 


depends only on er and 6. To the best of our knowledge, 
Hamiltonians of this type were previously studied for a 
single exciton [20] |BB], and very recently for a trion [4B], 
but not for a biexciton or any other bound state of more 
than three carriers. Denoting the lowest eigenvalue of 
hs(cr) in Eq. (13) as fs(cr ), we obtain the general result 


for an arbitrary number of carriers 

fs{o) 


£s ( cr,r 0 ) = - 


2 f 0 

Ne + N h — (N e - N h ) 2 

4fn 


- log 4fo -7 


(15) 

This expression suggests that by finding a single eigen¬ 
value of Hamiltonian (13) for specific 5 and a, one can 


recover the entire dependence of £s(cr, f 0 ) on r o in the 
limit of strong 2D screening |20l HBI [BB] . 

Equation (15) becomes very simple in the case of a 
single exciton 


£x (&, fo) = 


log 4r 0 - 2 fx(cr) - 2j 
4 fn 


(16) 


This expression has a “built in” applicability range in the 
sense that it becomes unphysically negative at fo < 1 if 
we assume fx{<x) ~ 1. 


IV. PATH INTEGRAL MONTE CARLO 

In this section we briefly outline the general ideas be¬ 
hind the Path Integral Monte Carlo (PIMC) methodol¬ 
ogy and specific details of its application to finding many- 
body eigenstates of Hamiltonian ([I]). Most of the static 


properties of a finite-temperature system defined by a 
Hamiltonian operator H can be obtained from a gener¬ 
ating functional (partition function) |67| 


Z = Tr [e~ pH ] , 


(17) 


where /3 = 1 /ksT, and Tr represents the quantum me¬ 
chanical trace. Functional derivatives of the so defined 
partition function with respect to H can then be used to 
evaluate most of the relevant observables of the system. 
For example, the charge density of the system can be ob¬ 
tained by introducing an auxiliary electrostatic potential 
to H and then finding the functional derivative of Z with 
respect to this potential. 

The quantum mechanical trace in Eq. © can be rep¬ 
resented by a discrete path integral 


N 


Z = / Y[dx r 


(x m+1 \e- TH \x r 




(18) 


m= 1 


where r = fi/N is an imaginary time step, and x m spec¬ 
ifies the coordinates of the entire system at the m th mo¬ 
ment of imaginary time. Periodic boundary conditions 
are assumed, x N+1 = x 1 . If the time step is small enough 
(i.e., approaching the limit of continuous path integral), 
then the matrix elements in Eq. (181 can be approxi¬ 


mated by the so called symmetrized Suzuki-Trotter de¬ 
composition formula as [43 IBS] 


(x m+1 \e- TH \x m ) = Cexp j_ ^ TOi 


[U(x m+1 ) + U(x m )] | + 0(t 3 ), 

(19) 


where x™ is a coordinate of the i th particle at the m th 
moment of the imaginary time, C is the normalization 
constant and U(x) is the potential energy part of H. The 
partition function can now be evaluated approximately 
as an approximate multi-dimensional integral representa¬ 
tion of the continuous path integral, Eqs. (18) and (19). 


Various observables of the system can be obtained 
by finding the functional derivative of the integrand in 
Eq. (18) with respect to H and then evaluating the re¬ 


sulting integral. One of the most efficient ways to eval¬ 
uate such integrals for large N = /3/r is to employ 
the Metropolis-Hastings Monte Carlo algorithm 
hence the name “Path Integral Monte Carlo”. 

It can be demonstrated that the error of the order • 


in Eq. (19) transforms into the error of the order of ~ t z 


for observables. Therefore, one can expect the following 
relation between the exact value of an observable A and 


its approximation An obtained using Eq. (19) 


B 2 

A »~ A(X ^2- 


( 20 ) 


If A is temperature-independent at very low tempera¬ 
tures, then this expression can be used to extrapolate 
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FIG. 1. PIMC results for the dimensionless exciton binding 
energy Ex as a function of p /N 2 are plotted as solid lines. 
System parameters are r o = 5 and <7 = 1. The extrapolation 
of the solid red and black lines to IV —> oo yields the binding 
energy of the exciton ground state, shown by dashed black 
line. Red and green horizontal dashed lines correspond to 
grid-based calculations with and without regularization of the 
inter-carrier interaction potential (see text for details). 


An to A by calculating An for a few values of N and /? 
and then performing a regression analysis. Figure[l]plots 
the dimensionless binding energy of a single exciton, Ex , 
as a function of f3 2 /N 2 , where the dimensionless inverse 
temperature is introduced as /3 = /3£u, where Eq is the 
Bohr energy introduced in Sec. [TT] As expected, the three 
solid lines corresponding to PIMC results for a range of N 
and p demonstrate a near-linear dependence when plot¬ 
ted against /3 2 /N 2 . Solid blue line corresponds to PIMC 
calculations with the largest temperature, /3 = 50. At 
this temperature, there is a noticeable admixture of ex¬ 
cited exciton states resulting in the apparent lower bind¬ 
ing energy (i.e., lower Ex)- The PIMC single exciton 
energies at lower temperatures, /3 = 100 and 150, are 
shown by the red and black lines, respectively. As can be 
seen, these two lines essentially coincide signifying that 
the temperature is low enough to sample only the ground 
state of the single exciton. Fitting these two lines with 
Eq. (201 produces the binding energy of a single exciton 
in the limit N —> oo, shown as the horizontal dashed 
black line. 

To benchmark the PIMC result, we have also per¬ 
formed a simple grid-based diagonalization of the single- 
exciton Hamiltonian. The result is shown by the hori¬ 
zontal dashed red line and is seen to be within the error 
bar of the PIMC result. 

It turns out that two modifications to Hamiltonian 
(|T|) are necessary to make it suitable for PIMC calcula¬ 
tions. The first modification concerns the fact that PIMC 
calculations are always performed at finite temperature 
(/3 < -poo), and at any finite temperature the entropic 
contribution to the partition function favors the complete 


dissociation of Coulomb-bound states. To prevent that 
we have included a weak attractive harmonic potential to 
each pair of charge carriers. Via the series of numerical 
tests we have demonstrated that at sufficiently low tem¬ 
perature it is always possible to tune the strength of this 
harmonic potential so that it is weak enough not to dis¬ 
turb the ground state properties of multi-carrier bound 
states, and strong enough to prevent dissociation. 

The second modification of Hamiltonian ([l]) concerns 
the singularity of the Coulomb interaction at vanish¬ 
ing inter-particle distances. Why this is a problem for 
PIMC calculations and how to solve it in the case of 
unscreened Coulomb interaction by introducing a reg¬ 
ularized Coulomb potential is discussed elsewhere 70] . 
Our generalization of the proposed solution consists of 
two steps. First, we perform a non-linear rescaling of 
the inter-particle distances when calculating the poten¬ 
tial energy 

Xij -p ) = [1 - exp(—apj/A)] -1 Xij. (21) 


The rescaled distance does not go below A, so the po¬ 
tential energy does not diverge. For the case of un¬ 
screened Coulomb interaction, this rescaling results in 
the regularized potential used in Ref. m- Second, one 
can notice that for potentials that are monotonic func¬ 
tions of inter-particle distances, e.g., Eq. ([2]), the in¬ 
troduced non-linear rescaling results in a non-vanishing 
perturbation of the total energy of the system already 
in the first order with respect to the difference be¬ 
tween potential energy functions with and without rescal¬ 
ing. This first-order deviation can be largely compen¬ 
sated by adding a well-behaved short-range potential 
to the regularized Coulomb potential, so that the inte¬ 
gral over the entire 2D plane of the difference between 
this resulting potential and the original potential van¬ 
ishes. Specifically, we substitute V(x) in Eq. |2]) with 
V'{x) = V(x') + Aexp(— x 2 /2X 2 ), where x' is defined by 
Eq. (21) and the prefactor A is found by requiring that 

J 0 fl xdx [V(x) — V'(x)\ approaches 0 as R —>■ oo. Since 
both x — x'{x) and exp(— x 2 /2X 2 ) decay exponentially 
fast at large distances, this integral approaches zero very 
quickly with R. This guarantees that the low-momentum 
spatial Fourier harmonics of V(x) and V'(x) are almost 
the same, so that V'(x) is a non-singular version of V(x) 
that is expected to reproduce the low-energy Coulomb 
scattering correctly. So introduced V’{x) is thus nothing 
but a local pseudopotential very similar to those used in 
atomistic quantum Monte Carlo methods m- The grid- 
based calculations with and without these two modifica¬ 
tions to the interaction potential (i.e., harmonic potential 
and zero-distance regularization) are shown in Fig. [l] by 
dashed red and green lines, respectively. As can be seen, 
the relative deviation of one from the other is ~ 3 x 10 -3 . 
In all the PIMC calculations in the present work we kept 
the deviation at this level or less by adjusting A. 

In general, the accurate treatment of the fermion 
statistics is required to correctly reproduce eigenstates of 
Hamiltonian |l]) when more than a single charge carrier 
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of the same type (charge) is present. However, appli¬ 
cations of quantum Monte Carlo methods to fermionic 
systems are plagued by the so called “fermion sign prob¬ 
lem” 071 EH Eg. The easiest way to avoid this prob¬ 
lem is to treat charge carriers as distinguishable parti¬ 
cles and not fermions. This could be done if carriers of 
the same type do possess “flavors”, i.e., quantum num¬ 
bers that are conserved by all the interactions in the sys¬ 
tem. Then, if two carriers have different values of at 
least one of such quantum numbers, they become effec¬ 
tively distinguishable. For example, if the only present 
flavor is the spin, then two holes (or two electrons) with 
opposite spin projections can be treated as distinguish¬ 
able particles within a positively (or negatively) charged 
trion 031 $J tSS] ■ In the case of ML-TMDC or stanene, 
the two two-valued flavors are the spin and the pseu¬ 
dospin., where the latter labels the two degenerate K and 
K' valleys in the Brillouin zone [74]. These two flavors 
are independently conserved only approximately since 
strong excitonic effects and valley-dependent spin-orbit 
coupling could in principle lead to mixing between the 
spin and pseudospin projections. However, this mixing 
was found to be weak due to the large momentum and 
high energy barrier separations of the two valleys m- 
At these conditions, the “flavor multiplicity” is four be¬ 
cause up to four same-charge carriers can be distinguish¬ 
able in ML-TMDC and stanene. Even more, same-charge 
carriers can become distinguishable in lead chalcogenide 
semiconductors where there are four degenerate valleys, 
and, therefore, the total flavor multiplicity (i.e., including 
spin) is eight [75M77] . 

Pauli repulsion suggests that the lowest-energy con¬ 
figuration of a multi-carrier bound state would consist 
of effectively distinguishable particles (i.e., with different 
spin or pseudo-spin projections), if the number of same- 
type particles is less than the overall flavor multiplicity, 
which is four in ML-TMDC and stanene, and eight in 
lead chalcogenides. Since, this work is focused on prop¬ 
erties of lowest-energy configurations of X, T+, T_ and 
XX, we neglect the fermion statistics and treat all charge 
carriers as distinguishable. 


V. NUMERICAL RESULTS 

PIMC-generated radial distribution functions for a sin¬ 
gle exciton, trion and biexciton at cr = 1 and f o = 5 are 
plotted Fig. 0] Radial distribution function (RDF) for 
two particles (1 and 2) is introduced as 

g(n 2 ) = / \\_dxi \'S/(x 1 ,x 2 , ...,x N )\ 2 S(r 12 - \xi - x 2 \), 

( 22 ) 

where ^(aq, x 2 , ...,Xn) is the wave function of the ground 
state of the Hamiltonian in Eq. ([I]). So introduced radial 
distribution function is normalized as f n dr 12 g{r 12) = 
1 , and since the wave function is finite everywhere, the 
radial distribution functions demonstrate the standard 



FIG. 2. RDFs for inter-carrier distances within an exci¬ 
ton (black), trion (red) and biexciton (blue) for a = 1 and 
fo = 5. Electron-hole and hole-hole RDFs are shown in solid 
and dashed lines, respectively. Distance (horizontal coordi¬ 
nate) is given in atomic units. RDFs are normalized so that 
f df g(r) = 1. 


2D behavior as they become linearly proportional to ri 2 
at small distances. 

Of all the RDFs in Fig. [2j the excitonic one (black 
line) is seen to result in the smallest mean electron-hole 
distance. This is of course expected as trions and biexci- 
tons are supposed to be much “looser” bound states with 
larger mean inter-carrier distances [621166] . If we assume 
that fo = 5 is large enough so that the strong 2D screen¬ 
ing limit is valid and Eq. (fTTI) is accurate, then we can es- 

■ ^ ~ _1 /2 

timate the excitonic size as (r x ) ~ \fx,eh) ~ ? 0 (ix, eh )- 
Furthermore, since all the parameters in h$ in Eq. (13) 
are on the order of ~ 1 at a = 1, we can estimate the 
exciton size in ^-coordinates as (£x,eh) ~ 1. This results 
in {fx) ~ 2, which is in a rather good agreement with the 
position of the maximum of the excitonic RDF in Fig. [2] 

Binding energies of a trion and biexciton are conven¬ 
tionally introduced as 


trb 

'T 


and 


£ xx = £ v x - 2£x, 


’■X, 

(23) 

2£x, 

(24) 


that is the binding energy of trion is defined as the energy 
released upon adding one more carrier to an exciton. The 
biexciton binding energy is defined as the energy released 
upon the formation of a bound state of two single initially 
isolated excitons. At f 0 —> 0, ££ and £ xx are typi- 
cally much smaller in magnitude than the single exciton 
binding energy £ x (£x = £ x in this work) [57J [5SJ 152] - 
IMj . Below we demonstrate that this qualitative picture 
remains true at arbitrary fo- As a result, carriers can 
be thought to form isolated excitons within the zeroth- 
order picture. Consequently, the first-order correction to 
this picture consists of accounting for weak interactions 
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FIG. 3. Dimensionless energy of an exciton, trion and biex- 
citon, £s((J,fo), as a function of fo at a = 1. Thin dashed 
black line represents the results of the grid-based calculations 
of the single exciton energy. Energy of a dissociated biexciton 
(i.e., two isolated excitons), obtained from grid-based calcula¬ 
tions, is shown by the thin dashed blue line. PIMC results for 
the exciton energy are shown by the thick dashed line. Red 
circles and blue squares represent the PIMC calculations for 
trion and biexciton energies, respectively. 


of these neutral bound electron-hole pairs between each 
other and with uncompensated charge carriers (e.g., in a 
trion) [58]. For instance, an electron in a positive trion 
is expected to “see” one of the two holes being “exci- 
tonic” with characteristic electron-hole distances being 
comparable to the size of a single isolated exciton. The 
other hole is farther away as it forms a much more weakly 
bound hole-exciton complex. These heuristic considera¬ 
tions are in agreement with what can be deduced from the 
electron-hole RDFs for a trion (solid red line) and biex¬ 
citon (solid blue line) in Fig. [2j Indeed, each of these two 
lines has two prominent features: (*) a maximum very 
close to the excitonic one (solid black line), thus repre¬ 
senting the “excitonic” hole, and ( ii ) a heavy shoulder 
at larger distances, representing the other more weakly 
bound hole. The Coulomb repulsion between carriers of 
the same charge results in the maxima of the hole-hole 
RDFs (dashed red and blue lines) being displaced toward 
larger distances, as compared to those of the electron-hole 
RDFs. 

Dimensionless energy of an exciton, trion and biexci¬ 
ton, £s(<j,ro), as a function of fo at a = 1 is shown in 
Fig- [3] The results of the grid-based calculations of the 
single exciton energy are given by a thin dashed black 
line and seen to be in excellent agreement with the cor¬ 
responding PIMC results (thick dashed black line). Red 
circles and blue squares represent PIMC results for trion 
and biexciton energies, respectively. As can be seen, the 
trion energies (red circles) are very close to the single ex¬ 
citon ones (thick dashed black line). Similarly, biexciton 
energies (blue squares) are very close to the energies of 


fo 

£x 

a x 

St 

CTt 

£xx 

<?xx 

0 

1 

- 

1.125 

- 

2.1928 

- 

0.015 

0.914 

2e-3 

1.024 

3e-3 

1.993 

5e-3 

0.05 

0.783 

2e-3 

0.873 

3e-3 

1.692 

6e-3 

0.5 

0.385 

le-3 

0.4184 

2e-4 

0.805 

le-3 

5 

0.10569 

2e-5 

0.1132 

8e-5 

0.2169 

6e-5 

50 

1.992e-2 

2e-5 

2.104e-2 

le-5 

4.049e-2 

2e-5 

500 

3.072e-3 

3e-6 

3.203e-3 

2e-6 

6.213e-3 

2e-6 


TABLE I. PIMC results for dimensionless energies of single 
excitons {£x), trions (£t) and biexcitons (£x x) calculated 
for a range of fo at a = 1. Error bars for these energies are 
given by ax, ar and axx, respectively. The first column 
is the dimensionless 2D screening length fo. At vanishing 
2D screening (fo —> 0), the trion and biexciton energies are 
obtained from Refs. [63] and [56], respectively. The analytical 
solution for the exciton binding energy at fo —> 0 is £x = 1 

m- 



fo 

FIG. 4. Binding energies for single excitons (black), trions 
(red) and biexcitons (blue) vs fo at a = 1. PIMC results are 
shown as solid lines with error bars. Large-fo asymptotics, 
Eq. ( |15| ), are shown by dashed lines. Red and blue filled 
circles represent the trion and biexciton binding energies at 
fo —^ 0 obtained from Refs. [56| [531 . 


two isolated excitons (thin blue dashed line). PIMC re¬ 
sults in Fig. [3] constitute the main result of this paper 
since, with the appropriate rescaling using Eq. ©, they 
can be straightforwardly used to evaluate the exciton, 
trion and biexciton energies for an arbitrary 2D semicon¬ 
ductor if the electron and the hole effective masses are 
isotropic and not too different from each other. To fa¬ 
cilitate the use of this data in analysis of experimental 
results or benchmarking of theoretical models, we also 
provide this same data in Table [T] 

Binding energies of trions, Eq. ( [23] ), and biexcitons, 
Eq. (24), are compared to the single exciton binding en¬ 
ergy, £x, in Fig. i Specifically, binding energies of sin¬ 
gle excitons, trions and biexcitons as a function of fo (at 
a = 1) are depicted by black, red and blue solid lines, re- 
































spectively. The lowest 2D screening length for which we 
have performed PIMC calculations is fo = 0.015. PIMC 
calculations with even smaller f 0 become increasingly dif¬ 
ficult because the Coulomb potential becomes more sin¬ 
gular at vanishing inter-particle distances EH- Fortu¬ 
nately, the numerical results for the trion and biexciton 
binding energy at exactly f 0 = 0 were obtained previ¬ 
ously by other methods [56, 63]. These energies are de¬ 
picted in Fig. [4]by red and blue filled circles, respectively. 
The analytical solution for the 2D exciton, Ex (cr, fo) = 1, 
is depicted by a black filled circle. As can be seen, our 
numerical results for a few smallest ro come quite close 
to the fo = 0 limit, so that a relatively smooth interpo¬ 
lation is possible and there is no need for explicit PIMC 
calculations with very small 2D screening lengths. 

At very large f 0 , Eq. (151 is expected to become accu¬ 
rate at describing energies of multi-carrier bound states 
in 2D semiconductor materials. We performed PIMC 
calculations for the ground states of hs, Eq. (131, for the 
exciton, trion and biexciton at a = 1. The numerical 
results are 


fx = 0.1800 ±0.0005, 
f T = 0.0382 ± 0.0005, 
fxx = 0.283 ± 0.002. (25) 


Up to numerical accuracy and notation, fx obtained here 
coincides with the one obtained recently in Ref. [46]. Us¬ 


ing these constants in Eq. (151, we plot the resulting 
asymptotic dependencies of binding energies as dashed 
lines in Fig. [4] As can be seen, a good agreement be¬ 
tween the exact PIMC results for the single exciton bind¬ 
ing energy and the asymptotic dependence, Eq. (161, is 
already reached at fo ~ 10, and becomes progressively 
better at larger fo- On the other hand, the asymptotic 
dependence deviates significantly from the exact solution 
at 1 < fj < 10, finally becoming unphysically negative 
at fo < 1. 

Evaluation of the trion and biexciton binding energies 
in the lim it of l arge f o is done via combining Eq. (15) with 
Eqs. (|23|) and (24). This produces £$, = (fx — Jr) /2fo 


and Ef x = (2 fx- fxx) /2f n , he., the /^-independent 
terms in Eq. (15) cancel each other out when Ef and 
E b xx are evaluated. Consequently, these expressions yield 
two straight red and blue dashed lines for trions and 
biexcitons, respectively, when plotted in the log-log axes 
in Fig. [4] The PIMC results for a trion and biexciton 
are seen to converge slower with fo to their correspond¬ 
ing large-fo asymptotics, than those for a single exciton. 
The rationale is that the inter-carrier distances within 
trions and biexcitons are larger than those within single 
excitons (see Fig. [ 2 ]), so r 0 has to be larger to provide 
the same extent of agreement between the finite-fo and 
asymptotic f 0 — > 00 results. 

In Fig. [4] the dependence of binding energies on fo is 
qualitatively the same decaying function for excitons, tri¬ 
ons and biexcitons except for a slowly varying prefactor 
seen as an almost constant vertical shift in the log-log 


-2 -1 0 1 2 3 

10 " 10 10 10 10 10 



FIG. 5. (a) PIMC results for the Haynes factor - binding 

energy of a trion (solid red line) or biexciton (solid blue line) 
normalized by the exciton binding energy - as a function of 
screening length fo at a = 1. Dashed lines of respective colors 
are the large-fo asymptotics, Eq. (151. Red and blue filled cir¬ 
cles represent the trion and biexciton Haynes factors at f 0 —¥ 0 
obtained in Refs. [5511531 . Red triangles are the trion Haynes 
factors, obtained in Ref. m for a series of ML-TMDC mate¬ 
rials. Insets show schematically the energy ordering of single 
exciton, trion and biexciton peaks in a photoluminescence 
spectrum for weak (left) and strong (right) 2D screening. 

(b) The ratio of biexciton and trion Haynes factors as a func¬ 
tion of screening length fo. The black dashed line is a guide 
for the eye for r/xx/riT = 1. The blue filled circle represents 
the ratio of Haynes factors at f 0 — > 0 1561 l63l . 


axes when going from single excitons to trions or biexci¬ 
tons. To focus on this prefactor it is convenient to intro¬ 
duce the so called trion (biexciton) Haynes factor [78] as 
the binding energy of a trion (biexciton) normalized by 
the binding energy of a single exciton taken at the same 

ro 

Vt = (Et — Ex) I Ex, (26) 

and 


Vxx = (Exx - 2E x )/Ex- 


(27) 


The Haynes factors for a trion (solid red line) and a biex¬ 
citon (solid blue line) are plotted in Fig. [5][a) at cr = 1. 
Dashed red and blue lines correspond to the large-f 0 
asymptotic behavior, Eq. (151, for trions and biexcitons, 
respectively. Red triangles represent trion Haynes factors 
for a series of four ML-TMDC materials, calculated in 
Ref. [23] by means of a simple variational ansatz approx¬ 
imation for a trion wave function. As expected, these 
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variational results underestimate the exact binding en¬ 
ergy (solid red line), but by no more than ~ 20 — 35 %. 
On the other hand, it is seen that in the range of f 0 cor¬ 
responding to the red triangles, the large-fo asymptotics 
(red dashed line) can overestimate the exact trion bind¬ 
ing energy (red solid line) by up to a factor of ~ 2. The 
dielectric screening by the environment is expected to de¬ 
crease fo even further [ 22 , thus potentially making the 
large-fo limit even less accurate when applied to realistic 
ML-TMDC. 

Numerical results shown in Figs.[4]and[5](a) can be used 
to obtain the relative spectral positions of photolumines¬ 
cence (PL) peaks pertaining to a single exciton, trion and 
biexciton. Indeed, if one assumes that the final state of a 
radiative recombination of a trion is a single charge car¬ 
rier occupying the lowest-energy band edge state, then 
the spectral energy gap between single exciton and trion 
PL peaks, Hcox — Hut, is identical to the trion binding 
energy [79]. Similarly, assuming that the final state of a 
biexciton radiative recombination is a single exciton in 
its lowest energy configuration, huix — Swx.y is identical 
to the biexciton binding energy. 

Figs. [4] and [5](a) suggest that the binding energy of a 
trion is larger than that of a biexciton at small fo, while 
the opposite is true at large fo- The crossover, i.e., where 
the trion and biexciton binding energies become identical, 
occurs at fo ~ 0.5. This crossover is further highlighted 
in Fig.gb), where the ratio of Haynes factors for a biex¬ 
citon and trion is plotted as a function of fo- Left and 
right insets in Fig. [5](a) show schematically the relative 
position of A', T and XX PL peaks at small and large fo, 
respectively. The former (left) PL spectrum is typical for 
semiconductor quantum wells where the effects of dielec¬ 
tric confinement are weak and the the dielectric screening 
of the Coulomb interaction is local (fo —> 0) [5T] ■ On the 
other hand, note that the biexciton binding energy can 
in principle be smaller than that of a trion, as has been 
demonstrated both experimentally and theoretically for 
semiconductor carbon nanotubes [5DH52 - In particular, 
a crossover very similar to that in Figs. [3] and [5] has been 
theoretically observed in Ref. [80] for single-wall semi¬ 
conductor carbon nanotubes. The crossover has been 
attributed to non-local screening effects, which is again 
very similar to what is observed in the present work. 


A. Energy variation with a 

All the calculations above were performed assuming 
identical effective masses of electrons and holes, i.e., 
(7 = m e /mh = 1. This is because the electron and 
hole effective masses are similar in ML-TMDC, a > 0.5 
mmm, which is presently one of the most important 
classes of 2DS. In what follows we show that the varia¬ 
tion of binding energies of multi-carrier bound states is 
typically not very strong at a ~ 1. Under these condi¬ 
tions, a binding energy at cr = 1 could be a reasonable 
approximation for binding energies at a > 0.5. 



a 


FIG. 6 . The dependence of the positive trion (solid red line), 
negative trion (dashed red line) and biexciton (solid blue line) 
energies on a at fo = 5. The energies £s(cr) are normalized 
by £s (cr = 1). Inset shows the same data, emphasizing the 
weak variation of £<s(<t) with cr at cr ~ 1 . 


The weak variation of biexciton energies with cr at 
a ~ 1 and arbitrary fo was already suggested in Sec. 13 
on the basis of vanishing derivative d£xx(c, ^o) / d<j\ a=1 . 
Single exciton energy Ex does not depend on the mass 
ratio whatsoever. That trion energies also vary weakly 
with a at cr ~ 1 is not known a priori, thus calling for 
direct numerical testing. Figure [6] shows the PIMC re¬ 
sults for cr-dependences of multi-carrier bound state en¬ 
ergies normalized to the respective energies at a = 1, i.e., 
£${&, fo)/£a(l, fo)j 5 = T + , T_ and XX. The 2D screen¬ 
ing length is set to f 0 = 5. As can be seen, the deviations 
of the biexciton and trion energies from those at cr = 1 do 
not exceed ~ 0.3% for cr > 0.5 (see the inset). Similarly 
weak variations of trion energies at cr ~ 1 were previously 
observed in the limit of the vanishing (fo -A 0) and very 
strong (f 0 -A 00 ) 2D screening [151155]- It is thus not 
unreasonable to expect a universally weak variation of 
Eg (a, fo) with cr at arbitrary f 0 if cr ~ 1. It has to be em¬ 
phasized, however, that the absolute energy of a multi¬ 
carrier bound state, Eq. does depend on effective 
masses not only through cr = m e /mh , but also through 
the reduced mass /i = (m“ 1 +m)[ ) -1 . Consequently, the 
absolute energies can vary strongly with masses even if 
cr remains the same. However, these variations are “triv¬ 
ial” in a sense that no new numerical analysis is needed 
since /.1 enters Eq. § algebraically. 

The dependence of £s(cr, f 0 ) on cr becomes stronger at 
cr < 0.5 for trions and biexcitons. This is especially no¬ 
ticeable for a positive trion and a biexciton, where two 
holes become very heavy at cr -a 0. Under these con¬ 
ditions, the adiabatic Born-Oppenheimer approximation 
can be applied, so that the first fmite-cr corrections to 
E,xx(cr = 0,fo) and £t + (& = 0,fo) can be calculated as 
a zero-point energy of vibrational motion of the two holes 
relative to each other, Eq. |To|). Resulting non-analytical 
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FIG. 7. Biexciton hole-hole RDF for a range of a at fo = 5. 


dependence, £,5(0, fo) — Es(<r, fo) oc ct 1 / 2 , is consistent 
with the rapid growth of the positive trion and biexci¬ 
ton energies when a approaches zero. 

The emergent “classicality” of heavy holes can also be 
seen in Fig. [7] where the hole-hole RDF becomes progres¬ 
sively narrower when the effective hole mass increases. In 
the limit of <7 — >■ 0, the two holes would become classical 
particles and, therefore, the hole-hole RDF would be a 
delta-function localized at finite f at zero temperature. 


VI. MATERIAL EXAMPLE: WS 2 

The numerical results obtained in this work can be 
directly used to analyze experimental data. For conve¬ 
nience, we compiled all the numerically calculated ener¬ 
gies of multi-carrier bound states in Table]]] Interpolating 
between the data points in this table and using Eq. [g] al¬ 
lows one to predict binding energies of excitons, trions 
and biexcitons. Other than data in the table, the follow¬ 
ing parameters of the system are required: n - the effec¬ 
tive dielectric constant of the environment, % 2 D - the 2D 
polarizability of 2DS, as well as m e and - the effective 
masses of charge carriers. The two-dimensional polariz¬ 
ability and effective masses can be obtained from the 
electronic structure theory calculations [52i 23] . Some or 
all of these parameters can also be obtained from exper¬ 
iment [55j . 

To provide a concrete example, we demonstrate here 
how binding energies of multi-carrier bound states could 
be obtained for a specific material at realistic condi¬ 
tions - a monolayer of WS2 put on top of a SiC>2 sub¬ 
strate. The 2D polarizability of WS 2 is Xid = 6.03 A 
[ 53 ] . The optical dielectric constant of the SiC>2 sub¬ 
strate is Kg? q 2 =2.1, resulting in the effective dielectric 
constant of re = («gfo, + K a i r )/2 = 1.55. Then, from 
Eq. ([3]) we obtain tq ~ 24 A. Using the exciton re¬ 
duced mass [i = 0.16 [ 53 ], we obtain the Bohr energy of 
Eq « 1.81 eV and the Bohr radius of a ss 0.51 nm, which 


results in fo ss 4.77. Assuming er = 1 and interpolating 
between the data points in Ta ble [T] produces Ex = 0.109, 
Et = 0.117 and Exx = 0.224 [54f l Using Eq. , these 
values are straightforwardly converted to exciton, trion 
and biexciton binding energies of 0.39 eV, 28 meV and 21 
meV, respectively. Inequality of electron and hole effec¬ 
tive masses in WS 2 , a ~ 0.7 — 0.85 [35]US], can be shown 
with the help of Fig. [6] to change these binding energies 
by at most ±2% of their respective values. As discussed 
in Sec. |III A] the exciton binding energy does not depend 
on a whatsoever, as long as the reduced mass remains 
fixed. The obtained exciton and trion binding energies in 
WS2 are consistent with the experimental values of 0.32 
eV [22] and ~20-40 meV, respectively [22] |3T] . Changing 
the effective dielectric constant of the environment can 
strongly affect these binding energies (see e.g., Ref. [35] L 
For example, substitution of SiC>2 with HfCU [86; changes 
the optical dielectric constant of the substrate from 2.1 
to «4 m- The dimensionless screening length then be¬ 
comes fo ~ 1.83, yielding the exciton, trion and biexciton 
binding energies as 0.27 eV, 21 meV and 18 meV, respec¬ 
tively. 

This same example reveals a peculiar contradiction be¬ 
tween numerical results obtained in this work and very 
recent experimental observations of biexcitons in ML- 
TMDC. Specifically, the dimensionless screening length 
for realistic ML-TMDC-based devices can be estimated 
to be r 0 « 5 [221 [36]. Under these conditions, accord¬ 
ing to Figs. [4 and [5] the ordering of PL peaks is ex¬ 
pected to be similar to that shown schematically in the 
r.h.s. inset of Fig. (5ja), i.e., the trion peak is expected 
to be lower in energy than the biexciton peak. The ex¬ 
perimental observations suggest otherwise [35] [38], i.e., 
the biexciton PL peak is located at lower energies than 
that of a trion. More quantitatively, the ratio of bind¬ 
ing energies of a biexciton and trion is estimated from 
our calculations to be r) X x/v t ~ 0.74, Fig. [5](b). On 
the other hand, the value of « 1.65 for the same ratio 
of binding energies can be extracted from Ref. [36j. This 
latter value can only be approached with our numeri¬ 
cal results if we assume a strictly local screening, i.e., 
fo —> 0. However, such a local screening cannot explain 
the non-hydrogenic Rydberg series for an exciton in ML- 
TMDC [55]. Provided that the experimental PL peak 
assignment and binding energies are correct, we need to 
critically rethink the applicability of the simplistic effec¬ 
tive mass model, Eqs. |T]) and ([2]), for describing multi¬ 
carrier bound states in ML-TMDC. Further research in 
this direction is thus required. 


VII. CONCLUSION 

In this work, using the Path Integral Monte Carlo 
methodology we have numerically analyzed the problem 
of multi-carrier bound states in 2D semiconductors where 
the electron and hole effective masses are isotropic and 
not very dissimilar. Specifically, we have calculated en- 
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ergies of single excitons, trions and biexcitons in their 
lowest-energy configurations for a range of dimensionless 
2D screening lengths, f o- These energies were demon¬ 
strated to converge correctly to the previously analyzed 
limits of vanishing 2D screening (trions and biexcitons) 
[56] [63], as well as very strong 2D screening (trions) [46] . 
The latter reference could actually be considered com¬ 
plementary to the present work, since with the help of 
Eq. ([6]) the entire problem can be reformulated in terms 
of an unknown function of two independent parameters: 
the mass ratio, a = We/m^, and the dimensionless 2D 
screening parameter, f o- In the present work we have sys¬ 
tematically investigated the dependence of this function 
on the second parameter assuming a ~ 1. The depen¬ 
dence on the first parameter was investigated in Ref. [46j 
in the limit of very large fo- 

Two limitations of the PIMC methodology, as imple¬ 
mented in the present work, are (i) the effective mass 
approximation, and (ii) its inability to accurately access 
energies of excited states of multi-carrier bound states. 
The effective mass approximation is typically justified if 
an exciton is of the Wannier-Mott type, i.e., the average 
distance between the charge carriers is much larger than 
the size of the unit cell of the semiconductor material. 
As discussed in Introduction, this is satisfied for ML- 
TMDC. In general, however, if an exciton is so strongly 
bound that the effective mass approximation fails, Bethe- 
Salpeter-like methods have to be used to assess excitonic 
effects [2U EH H2 138] 139] HQ. 

As for the second limitation, only excited states of sin¬ 
gle excitons have been experimentally observed to date 
mma, and it is straightforward to obtain excited sin¬ 
gle exciton states via grid-based numerical methods (see 
this work and Ref. [23]). If excited trions and biexcitons 
have to be analyzed, one can use density matrix-based 
methods [42 J [45] or a generalized PIMC [88] . 

Finally, we would like to emphasize that even though 
ML-TMDC constitutes one of the most promising classes 
of 2D semiconductors at present, the results of this 


work are very general and could be applicable to many 
other materials and systems as well. For example, lead 
chalcogenide nanosheets are effectively 2D semiconduc¬ 
tors with almost identical isotropic electron and hole ef¬ 
fective masses (39) [99] ■ These materials do also possess a 
large “flavor multiplicity” of eight (see Sec. IVI required 
to avoid the fermion sign problem when applying PIMC 
to multi-carrier bound states in their lowest energy con¬ 
figurations. Other potential applications include stanene 
nmiH] and ultrathin organic-inorganic perovskite crys¬ 
tals [T9]. In particular, the latter materials have been 
shown to have large exciton binding energies |l9l com¬ 
plemented by a non-hydrogenic exciton Rydberg series 
suggesting a non-local screening of the Coulomb interac¬ 
tion, Eq. ©• Such perovskite crystals could have not too 
different electron and hole effective masses mm , ren¬ 
dering the numerical results of the present work directly 
applicable to them. 

Numerical results obtained here are not applicable to 
phosphorene 631 - 66] due to a very large anisotropy of ef¬ 
fective masses. However, the PIMC methodology can be 
straightforwardly modified to account for this. Analysis 
of the multi-carrier bound states in 2DS with the large 
anisotropy of effective masses will be the focus of our 
future work. 
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